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The random motion of a particle with nonlinear damping is investigated. 
The spectrum oj the velocity of the particle is obtained by solving the as- 
sociated n on stationary Fokker-Planck equation and also by using the 
equivalent-linearization technique. The first procedure yields an exact 
solution in terms of Laguerre polynomials. The second leads to simple, 
approximate results which are valid for cases where the small nonlinearity 
assumption holds. Results obtained by these two methods are compared 
and good agreement is observed over a large frequency range. 

I. INTRODUCTION 

The recent advance in space and communicational technologies has 
led engineers to numerous difficult but fascinating problems in regard 
to the structural dynamics in random environments. For example, 
Hempstead and Lax have investigated noise in self-sustained oscilla- 
tion j 1 ' - Ariaratnam and Sanker have studied the dynamic snap- 
through of shallow, arch-type aircraft components under stochastic 
pressure. 3 In this paper the random vibration of a simple mass with 
nonlinear clamping is studied. The nonlinearity of the system is in- 
troduced to linear viscous damping by adding to it an extra term 
which is inversely proportional to the first power of the current ve- 
locity. Emphasis of the analysis is placed on finding the power spectral 
density of the random motion. 

Two different approaches are used to obtain the desired solution. 
First, the exact spectrum is found by solving the associated nonsta- 
tionary Fokker-Planck equation in terms of the eigenfunction expan- 
sion of the degenerate ordinary differential equation. Second, approxi- 
mate solutions are obtained by using the equivalent linearization 
technique by which the original nonlinear system is converted to an 
equivalent linear one. The equivalent linear system, constructed by 
the least mean square error criterion and based on the small non- 
linearity assumption, is then solved by standard linear theory. 
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II. NONSTATIONARY FOKKER-PLANCK EQUATION 

Consider the first-order nonlinear system described by the follow- 
ing differential equation : 

x + F(x) = f(t), 0) 

which may be thought of as the velocity equation of a unit mass with 
nonlinear damping F (x) subject to a force f(t). 

Let us discuss the problem of obtaining the power spectral density 
of x(t) when f(t) is a random process. We limit the discussion to 
the case where f(t) is a stationary white gaussian process with the 
first two moments defined as 

</(fl) - (2) 

and 

(f(ti)Kt 2 )) = 2tts «(r, - t 2 ) (3) 

where s is a constant, the symbol ( ) indicates the ensemble average 
and 8 indicates the Dirac delta function. 

Caughey and Dienes 4 have investigated a similar problem for 
F{x) = k sgn x. We shall however consider a different case in which 

F(x) =Px-Z (4) 

< X < co 

where (3 is a constant and y is a smaller nonlinear coefficient. In case 
y = 0, equation 1 becomes the familiar linear differential equation. 
The Fokker-Planck equation which governs the transition proba- 
bility p{x | x, t) with given initial velocity x = x{t ) for the velocity 
x{t) at time t is 



v = Tx 



** "?>]+-*§ 



L 2 ? 

xfj ' " a °dx* 

where t = t — t . The initial and boundary conditions for equation 
5 are 



Vto = 5(x - x ) (6) 

and 

P(0, t) =p(«, = 0, (7) 
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respectively. As the time of passage t becomes sufficiently large, 
p(x | x, t) in equation 5 approaches a stationary value p s t(x) inde- 
pendent of t and initial condition 6. Setting p = in 5, such a station- 
ary density can be found by solving the degenerate stationaiy equa- 
tion. The result is: 

p. t = C exp \_-— ( 2 "" 7 log x ) (8) 

where C is the normalization factor determined by 



f 



p.t dx = 1. (9) 



The power spectrum is the Fourier transform of the autocorrelation 
function which is determined by the joint probability density p(x , 
x, t) . Thus from the relation 

p(x ,x,t) = p.,(xo)p{x | x, t), (10) 

we need to find the transition probability density, that is, the non- 
stationary solution of equation 5. Let p(x \ x, t) = T(t)X(x) in 5. 
It follows that 

T(t) + \T(t) = (11) 

and 

°>0 + i (sX)-!*!- (-)+£* -° (12) 

dx dx ts dx \x I j8 ' 

where a 2 = ir8 /p. If X,„ (x) is the eigenfunction and A m the correspond- 
ing eigenvalues satisfying equation 12 and the prescribed boundary 
condition 7, it can be shown that 8 

Pfel»,r)- ± X ~WffJ ,->-->. (13) 

m-0 P$t\Xo) 

In deriving 13, the following orthogonality condition has been used: 

Jxmzm j^-fc.. d4) 

Following the transformations adopted by Stratonovich, let /x = *4 
iy/vSo ~ 1), 2 = x 2 /2a 2 , and u = zr^X. Equation 12 becomes 

d 3 u du [ (1/2) + u + (X/2/3) (1/4) - ti' 1 

"a? + & + L * + — ? — J" = °- (15) 
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Equation 15 is a degenerate hypergcomctric differential equation which 
has eigenfunctions U„(z) with corresponding eigenvalues X„ = 2nfi: 

U% {z) = f* tvn 9-LF*to (16) 

where 

L;(,)=^e-z-"§,(e-^') (17) 

is the Laguerre polynomial of degree n. 

Transforming back to the original variables and applying equation 
14, we obtain the following normalized eigenfunctions: 

An{X) a [n! r(n + 2y. + l)r(2 M + 1)]* k ' 

From equations 13 and 10 the transition density and jointly density, 
respectively, can be found as 

, | N (.*) 2 C y> -k„ t,Z;-^n IS.J -2»g|T| / 1Q N 

pfo ' '■ r) = V "E™feT £«!r(« + 2„ + i) e (19) 

and 

p(x ,x,t) = —_ (zz n ) 



2)i+(l/2) -<2 + Io) 




a ti Ui WF*<ri r -wM 


(20) 


' n V n!r(n + 2 M + l)r(2 M + 1) 


2 2 

2 = ^ and z, = |V 





where 



From the above the autocorrelation function iMi-) of .x(0, where 

T = t — t ,\s 

rtoo /ieo 

■Rx(t) = / / p(x , x, t)x x dx dx 

°4 »^ + » + »! ' M (21) 
where 

/(,,,„) = f I" (zT + " m e-\z Y> + " /2) e-'°Ll 2 >\z)Ly{z )dzdz . (22) 
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In the appendix we show that 

f g2 „ + (i/2, -.- cm)/-, rf , _ ~r[2 M + (3/2)]r[n- (1/2)] 
.'o ' K {Z) dZ ~ 2(7r) J n! ' {2S) 

from which it follows that 

/( „, n) _ _L | r[2, + (3/2)]r[n - (1/2)] ^ (24) 

Substituting 24 into 21, we obtain the following expression for the 
autocorrelation function R x (t) : 

RM _L r 2 [2 M + (3/2)] f, r 2 [n - (l/2)]<f 2 " g|T| 

^ rJ 2tt T(2 M + 1) h (n!) 3 r(n + 2 M + 1) ' W 

which is a monotonically decreasing function of t. 

The power spectral density of x{t) can be derived from equation 
25 as the following: 

i M* r 2 [2 M + (3/2)] f, T 2 [n - (1/2)] 4n/3 

T W T(2 M + 1) £ (n!) 3 r(n + 2 M + 1) 4n 2 /3 2 + a. 2 " W 

Notice that S. T (o>) is again a monotonically decreasing function of 
w and has a spike at o> = 0. 

The nonstationary mean value of x{t) is given by 

(.f(r)> = P xp(x e \x,r)dx. (27) 

•'o 

Using equations 19 and 23 it can easily be shown that 

_ °T[2 M + (3/2)] f, T[n - q/2)]L™(z ) 2naiTi 
( ' l(T)) " - V2^ L.-0O » (n!) a r(n + 2 M + 1) e • (28) 

Because x = (2o- 2 2 ) % , we notice that (x(t)) depends on the initial 
velocity x . As t -» oo, (.t(t)) in equation 28 approaches its stationary 
value (.r) s , , which is given by 

= (2)'r[2 M + (3/2)] 

w " r(2 M + 1) *" (29) 

(x) st is independent of the initial condition .r . This stationary mean 
velocity, (x) 8 t, can also be found by using the stationaiy density 
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Pst(x) as follows: 



<x>.« = / xp.,{x) dx 

Jo 

= I xX {x) dx 

Ja 



(2)*<r 



r(2 M + l) 



)J> 



2M + (l/2) e -. L 2, (2) ^ 



_ (2)>r[2/z + (3/2)] 

r(2/i + i) 

The variation of <x(t)> is illustrated in Fig. 1. It has a maximum value 
x at r = and decreases exponentially to the stationary value (x)„ 
as given in equation 29. The nonstationary mean square value of x(t), 
given its initial value x 2 , is difficult to evaluate explicitly, but its 
stationary value <x 2 >.« can be found by integrating £«(«) in equation 
25 over the entire frequency range from — oo to + <» . By this procedure, 

T 2 [n - (1/2)] 
!) 3 r(n + 2 M + 1) 

(30) 

As expected, (x 2 ) st is also independent of the initial condition x . The 
variations of (x 2 ) at with the nonlinear coefficient k = y/irs are shown 
in Fig. 2. 



/ »N r Q t\ A "' r 2 [2 M + (3 /2)] f, _ 

(x ). t = J_ x BM rfo, = - r(2/i + 1} 2. (nl 






ti <x> st 



TIME, T 



Fig. 1 — Nonstationary mean velocity of nonlinear system subjected to random 
force. 



RANDOM MOTION VELOCITY 



2265 



S 



A<* 2 >t, 



EQUIVALENT LINEARIZATION 
USING TRUE DENSITY 




EXACT MEAN 
SQUARE VELOCITY 



^<X2> G! 



EQUIVALENT LINEARIZATION 
USING GAUSSIAN DENSITY 



0.02 



0.03 0.C4 0.05 0.06 0.07 
NONLINEAR COEFFICIENT, k 



0.09 0.10 



Fig. 2 — Comparison of mean square velocity nonlinear system subjected to 
random noise. 

in. EQUIVALENT LINEARIZATION TECHNIQUE 

The eigenfunction expansion of the degenerate ordinary differential 
equation of the governing Fokker-Planck equation is often difficult 
to find. In such cases it is sometimes convenient to use the perturba- 
tion 8,7 or the equivalent linearization 8,9 techniques to obtain desired 
quantitative results. If the nonlinearity of the system is small (that is, 
y/P « 1), these techniques provide the simplest means for obtaining 
approximate results. In the following, an equivalent linearization pro- 
cedure is used to derive the power spectral densities of the nonlinear 
velocity x(t) in equations 1 through 4. Let 



x + (3 t x + e(x) = f(t) 



(31) 



be the equation of motion of a system equivalent to that described by 
equations 1 and 4 in which /3 C is the equivalent linear stiffness and 



e(x) = fix + L - fi t x 
x 



(32) 



is the error function. If € (.t) is small and may be ignored, equation 
31 becomes a linear differential equation, and its spectrum can be 
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solved by standard technique. The equivalent linear stiffness p„ can 
be determined by the criterion that the mean square error is mini- 
mized. The mean square error is 

< € (z) 2 ) = J" 1(0 - fc)V + ^ + 208 - 0.)7jp(aO dx. 
Setting d(e(x) 2 )/d{3 c = we obtain 

/ $x 2 p(x) dx = / P.x*p(x) dx - \ yp(x) dx; 

J Jo J o 

therefore 

/3 / x a p(x) dx + y 



/ x a p(x) 

Jq 



dx 



-, + fr (33) 

The following two cases are considered: 
(i) Assuming p (x) is gaussian, that is, 

, . (2*\Y> ( (3x 2 \ 
p(x) = [- T ) exp^-^-J, 

then 

{x 2 )o = *l = =*■ (34) 

Substitution into equation 33 yields 

0..o = P0- +7/«.). ( 35 ) 

(ii) Using true (stationary) distribution, p gt {x) given by equations 
8 and 9 becomes: 

p., Or) = j w 



where 



= i/2 g= y ii+(T/ "- )i r r i + (7Ag.) "[ 



RANDOM .MOTION VELOCITY 2267 

Therefore the mean square value of x{t), using true distribution 
equation 36, is 

(x 2 ) r = I x 2 p.,(x) dx 

Jo 

= |[/; ( .)™exp(-^^)^ 

= ^(xs + 7 ). (37) 

Substitution into equation 33 yields 

P..r = «1 + o \ ) - fll + -*- - -£ + •■•)• (38) 

Comparison of equation 38 with equation 35 indicates that /? Ci o is 
the first-order approximation of /3 C , t- 

Now let us consider the simple linear system 

x + r x = f(t) (39) 

whose transfer function is given as 

According to the familiar linear theory of random processes, the 
power spectrum of x(t) in equation 39 is given by 

5 » = w+r^ (4i) 

where s is defined in equation 3. 

Substituting (3 C , G as given in equation 35 into equation 41, we obtain 

o 

Substituting /?,., T as given in equation 38 into equation 41, we obtain 

S ' Aa) = [^( 7 r 5o + 27)/(X + 7)] 2 +o J 2 - (43) 

By setting y = 0, both equations 42 and 43 give 



&(«) = sr 



+ « 2 
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which is the spectrum of the corresponding linear system. The mean 
square value and the power spectral density of the velocity response 
x{t), obtained by solving the Fokker-Planck equation, and the equiv- 
alent linear equations using gaussian and true distributions, respec- 
tively, are summarized in Table I. In this table k = y/vs is the non- 
linear coefficient. 

The mean square velocities (x 2 ), (x 2 ) , and {x 2 ) T are compared in 
Fig. 2. It is seen that for k < 0.035, that is, in a very small nonlinearity 
range, both linearization cases give higher mean square velocities than 
the exact solutions. For k > 0.035, equivalent linearization methods 
give larger results when using true distribution and smaller results when 
using gaussian distribution than the exact solutions. 

The power spectral density functions S x , a (u>) and S t , T (w) obtained 
by the equivalent linearization procedure, using gaussian distribution 
and true distribution of x(t), respectively, are compared in Fig. 3 in 
which B, = (P 2 /s )S x , a (o>) and B a = (0 2 /s o )S x . T (o,). Notice that both 
S x , a (u) and 8 x , T ((a) are monotonically decreasing functions, and the 
differences between them are negligible for small nonlinearity. 

The exact power spectral densities are compared with the equivalent 
linearized solutions in Fig. 4. For the exact solution of S x (u) as given 
by equation 25, the spike at co = is evaluated by normalizing 5,(w) 
to an area equal to that given by Fig. 2, that is, 

< x 2 ) = r 8 x (p>) du = 3.053 for k = 0.01. 

Jo 

Curves shown in Fig. 4 are seen to be monotonically decreasing. The 
equivalent linearized systems have higher power spectra of x{t) in the 
low- frequency region and lower power spectra of x(t) in the high- 
frequency region than the actual nonlinear system has. 

rv. CONCLUSION 

It has been shown that the exact expression for the two-dimensional 
nonstationary probability distribution of a class of simple nonlinear 
systems can be found in terms of the spatial eigenfunction expansion 
of the governing Fokker-Planck equation. Equivalent linearization 
techniques can be very useful in generating approximate response 
statistics for certain systems having small nonlinearities. In Figs. 2, 
3, and 4 good agreement has been achieved in the comparison of the 
exact and approximate mean square values and of the power spectral 
densities of the nonlinear random response. 
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Fig. 3 — Power spectral density by equivalent linearization. 
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Fig. 4 — Comparison of the exact and approximate power spectra. 
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APPENDIX 



Derivation of Equation S3 

The following formulae are used in the derivation of equation 23: 

f c~ at l p L:(t)dt 

(Re/3 > -1, Res > 0) (44) 
where 

F(a, 13; y; •) = 2 b\ (a, ft y; •) (45) 

is a generalized hypergeometric series which is defined as 

in which 

r(« + k) 



(46) 



WK '" 4 ' r( 7 -a)r( T -p) 

(Re 7 > Re (a 4- 0), Re 7 > Re p > 0). (47) 

Using equations 44 through 47, the integral involved in equation 
22 can be evaluated as follows : 









W* 


I» 


and 






(p)» 


r(p + *) 
r(p) 


A special case 


for 


equation 


45 is 


when z = 1 : 


bv~ /?•„,. 1 


N 


r( 7 )r( 7 


— a 


-p) 



r ^ (i/ v'li w (2) & 
•'0 

r[2 M + (1/2) + l]r(2 M + n + 1) 



n!r ( 2,4-l) F(-n,2 M 4-f;2 M + l;l) 
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_ r[2 M + (3/2)]r(2/i + n + 1) r(2 M + l)r[2/x + 1 + n - 2n - (3/2)] 
n! r(2 M + 1) r(2/x + 1 + n)r[2 M + 1 - 2 M - (3/2)] 

r(2 M + (3/2)] r [ n - (1/2)] = -r[2fi + (3/2)] Y[n - (1/2)] 
n!T[-(l/2)] 2(7r) J n! 

which is equation 23. 
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